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ABSTRACT 


tinea =e rnNOdSs 1TOregeneratring outcomes on multivariate 
normal random vectors are presented. A comparison is made 
memoctermine which method requires the least computer 
soul On time and memory Space when utilizing the 
Mermsco0/67. All methods wse as a basis a standard Gaussian 
random number generator. Results of the comparison study 
mmcai1ecacve that the method based on triangular factorization 
of the covariance matrix generally requires less memory 


merece and compuver time than the other two methods. 
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Pe CRODUCT EON 


Prequentily ls 1s desired to generate a sample of vectors 
derived from a specified multivariate normal distribution. 
Memecxample, this need arises in simulation studies in which 
Peemodel contains stochastic variates that are distributed 
mecordging to a multivariate normal. Since a computer is 
mascuently used to produce the sample vectors, the memory 
Space and computer time requirements become important cost 
factors when considering various methods for "generating" 
fmamdom vec ors. 

The problem of "generating" random vectors was discussed 
in 1948 by Herman Wold [Ref. 1]. Wold describes a method 
for construction of samples from a multivariate normal 
distribution. Wold used a method Bbascdwon triangulariZaurom 
of the covariance matrix. In 1962, Ernest M. Scheuer and 
Havid 5S. Stoller wrote a paper on the generation of multi- 
variate normal random vectors [Ref. 2]. They considered a 
method based on matrix equations and a method based on con- 
Gitional distributions. The first method uses a technique 
mmr lar to that of Wold's. Their conclusion was that the 
mathematics involved in the first method was simpler than 
that of the second method so the generation of random 
meceeors is Dest accomplished by the matrix equation 
eechnique. 

The objective of this paper is to examine the above two 


methods, hereafter called "triangularization" and 


TeodurLronineg respectively, and to discuss a third method 
which is referred to as the "rotation" method. Computer 
programs for each method are presented and a comparison of 
memory space requirements and execution times on the 
IBM 360/67 computer is made for the three methods. The 
format for making this comparison study is similar to that 
feca bY M. E. Muller in his paper on univariate normal 
generators [Ref. 3]. 

The following notation is used in this study: xX denotes 
Pepecuadom Variable whereas x is a real variable. A vector 
ts denoted by a bar under the letter. Listed below are two 


examples of the above notation: 


a al 
K X 
ne x= | ° 
xX ye 
Dp Dp 


Mites ols cuccuen Of the three methods, .eenerators for 
N(O,£) distributions are considered, where © is a general 
(positive definite symmetric) matrix. There is no loss of 
Benerality in assuming the mean vector to be the zero 
Beetor since a random vector Y distributed N(u,2) may be 
mem@erated by first generating X distributed N(0,Z) and then 


meking Y = X + u. 


Ii. METHODS 


Pee ROTATION METHOD 
Given a positive definite symmetric p x p matrix 2, the 


function 


-l] 
=o She eel 
f(x:5) = (2n)72P |r|7? eek 2 (1) 


is a p-variate normal density. A random vector X having 
fers distribution has mean 0 and variance covariance matrix 
y (Ref. 4]. Denote the distribution law corresponding to 
me density funetion (1) by N(0,2). If X (with p components) 
feedistributed N(Q,y), then Y = CX is distributed N(0,CrCc'), 
if C is nonsingular [Ref. 4]. The development of the 
febacvion method is based on this fact. If ap x p matrix C 
#s found such that Crc' yields a diagonal matrix, then 
elements on the diagonal of CrC' are the variances of the 
Berresponding components of Y. Since the covariance values 
moff diagonal elements of CEXC') are zero, the components of 
fare independent. 

For a symmetric matrix 2, there exists an orthogonal 
Meporix C such that C2C' = D where D contains the eigen 
Memues of 2 on its diagonal [Ref. 5}. Associated with the 
eigen values are eigen vectors which make up the rows of 
mee Matrix C. Since Y = CX is N(0,CzZC'), and there exists 
eee sich that ChC' is the diagonal matrix D, it follows 
mat the components of Y are independent normal variates 


meee Zero mean and variance values as given by the 


eorresponding eigen values of 2%. These variates may be 
generated using a univariate normal random number generator. 
Pemetnod Of generating the random vector Y is thus easily 
maeaolished, To obtain a fpenerator for the random vector 


X, a transformation is made using X = et 


eo Cee Cas 
manorthogonal matrix, this simplifies to the equation 
R= C'tY. Standard matrix multiplication of the orthogonal 


matrix C' and sample vectors generated on the random 


mecv0Or ¥Y yields the desired random samples of the vector x. 


ee CONDITIONING METHOD 
Penile var hewmen probability density function 1 may be 


written as: 


F(x) 546+ 9%) = £(x,|X55--- 5X) F(X5|Xgo+++ Xp) 


Cop one) (2) 


where P(x, |x5---) denotes a conditional density function. 


ft xX can be generated from f(x), then xX may be gener- 


-] 
aved by conditioning on Xo Xo may then be generated by 


Senditioning on Xn and X In @ similar manner, all the 


n-1° 


Pemponents of X may be generated. In what follows, this 
memeral approach is applied to multivariate normal 
mrs tributions. 


Clow unabevnenvyeevor A 1s partitioned into two sub- 


fectors x0) and 62) 


ager VecamernemecOvarlance Matrix Of X 


y and 2% then 


lee 2s Dias 
meecan be shown [Ref 4] that the conditional distribution 


ms correspondingly partioned into 2% 
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of oo given Ko? 7 x6?) is normal with mean 2,, aaa x62) 
-l 


and covariance matrix May - bio doo Boys Several defini- 


iemons are introduced atv this time as they are required in 


what follows. 


Let 
Ck k Ok kt] ie as 
"K+1,k tale ee ita ry 
pik). 
nk on k+l ——* enn 


and let Ue denote the cofactor of o,, in Sc) me taraeaie 
5 Sevoves Une scolacuor. Of aay in 2. Consider now the mean 
mame vertance Of the distribution determined by 

BO 1X y+ X,)- Suppose the covariance matrix t is 


partioned as follows: 


Pa Weel Os he reat renee Fr on! 
eee 
| 
eee ye ieee eee SOS ta a 
| 
se Nets el al aia) irate) = acai we) ge many EA 
DG a : ope 
| 
| 
| 
On k : Sterna On K+ ee OD 


(3) 


a 


The variance of f(x, [Xpypoee +X) is given by: 


Met pe SF 


V(X [Xp pors eX) = 2a, 7 Fo Foo Foy C4) 


To simplify this expression, the following theorem is used 
(Ref. 4]: |Z] = JZ,,. - 2,4 £55- 25,|-|2,,]. The square 

; ; iil ae 22 lS Oia 
matrix B55 is nonsingular since the covariance matrix 2% is 


positive definite. The variance (4) can be written as: 


5 (k) 
ME Nyy stem? FL) (5) 
The mean of f(x, | X47 5++ 9X) is given by: 
_ aie ey 
E(X | Xypgqere+ oy) = Zo Fon X 


Bere ct Substitution of components of the partitioned covari- 


ance matrix (3) into the latter expression yields: 


-1 
eS) (k+1), 


BCX Xyp pero %y) = Cy a eee ot nl © 


dare 


ene 


mesimplification of this expression can be made as follows. 


y(kt1) 7h 


First, examine ( By using standar.i matrix 


mrsepra, it can be shown that: 


2 


-1 
cp (kt) 77 


y (kt1) 
Ie el Sr cea all 
are, 
i KoPeee 
pp Okt) | 
(-1 )ktltny (k+1) 


Wear 


(k+1 ) 


Thy eT 


Ox41) 


ee ee 


yk+24nys (k+1) 


(-1 k+2 ,n 


Meme, pertorm the matrix multiplication 


(Oy tl «RK, k+2 


This yields a vector of dimension 1 x (n - k + 1). 


e] 


k,n 


-1 
; pp (kt1), 


+k+1. (k+1) 
(1 ene 


(-1 y ntk+2 »(k+1) 


n,kt2 


(-1) Pee)? 


The 


beanspose of this vector is the product of the factor 


1 


por, me 
rey (kt1) 0 v= } 


k+1,kt1°k k+l 


oe) 


be tk k+l 


[(~1) nykt+1°k k+l 


(k+1) 
Pea sie De ete 


(k+1) 
k+2 ,k+2-k,k+2 


k+lin, (k+1) (1 \Ktetn 


0 


re] 


(aE) 


an kt 2 


3 


k+ 1+n 


en iD) ea 


Dar ease 


al 


+2¢ni(kt]) 
(-1y" Diy nm Kn 


én, (k+1) 
n,n 


ie) 


Hie eel Bs 


(6 ) 


Mudcviply the vector (6) by [x 


yeqo+++o%,)] and the following 


results are obtained: 


E(X, |X) yq0°- 


oF (a2)x)] 


1 
Xn) So CerTY) 1882) Meta * (82) Mera 


where (a, ) indicates the first element of the vector (6), 


(a.) indicates the second element, etc. 


ion, it is 
In the same 


(k) 
Beek: 


as. 


The 


E(X, |x 


The use 


of es now 


- O 
mk tk 


dent normal 


Using expressions (5) and (7) for 9 


Kept? 


Upon close anspec— 


(k) 


apparent that (a, ) is equivalent to -2 


k+l ,k’ 
: (k) ; 

manner, (a. ) is a3) Tee ., and (a3) is 

mean of f (OG ecu can thus be represented 


) 1 (k) 


_ (k) (k) 
TT tRFTY) Pict Meta * Pet, ken! 


kt2,k*kt2 °**"°n,k*n 


(7) 


of these expressions in generating the components 


X 
ales 


described. Consider the. transformation 


My k = 


variates with zero mean and 


1,2,... n, where the Y's are indepen- 


unit variance. 


and HU 


k k respectively, 


foes transformation may be given explicitly as: 


1, 
meee) *, a a ye) 
k Ib Rell) k py (KFT) | Meade tale Ke Kiko Oe 


for k = 


ee 


(8) 


(k) 
nk *n 


co) 


Use of the equations (8) and (9) there- 


ici 


get me 


fore provides a method of generating samples of a random 


mecvor A, using univariate generators 


for the Y's. 
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C. TRIANGULARIZATION METHOD 
The covariance matrix, designated as 2, is a symmetric 
positive definite matrix. It can be written as 2 = TT! 


where T is a lower triangular matrix [Ref. 6], 


O11 %9 Cin on Ot i840 boy oil 
G5) % 0 eri ei wees ee 8 0 too - Uo 
onl On2 Cnn al Ae Oa 0 0 aa 


X 1s obtained by making the non-singular transformation 

X = TY where Y is distributed N(0,I). By previous remarks, 
KX is distributed N(0,TT'). The matrix T can be obtained 
from % by using the so called “square root method" [Ref. 6]. 


Bovuations for eo ingcermns of O,, are as follows: 


J 
en 0G)! Ja eee a ies A 
: rt 
i-l 5 

ee OVO Oh a ol Seed 

siveat ipod p=1 iP 
c =f cer Garey er (Rares 

Pe, = See Oa ; : a diam @! 

LJ O45 1J p=) %ip JP 
eee = 2.0) 8. Cpe SSS Aes ial 


1 
the vector Y consists of independent normal components with 
zero mean and unit variance. These may be generated using 
@ univariate normal random number generator. Standard 


Serix multiplication oi T and samples of Y yields samples 


Sie the random vector X. 


aS 


het ErnOD EVALUATION 


Computer runs were obtained using each of the three 
Mmeunods, with various size covariance matrices. Data con- 
Serming these runs is compiled in Tables I, I1, and IIL. 

A brief description of the procedure employed to obtain 
execution time and memory space requirements for the tables 
follows. The 2 tre programs can be subdivided into three 
meste portions. The first portion consists of the data 
teat must be entered into the computer. The second portion 
consists of all required steps, prior to the use of the 
basic univariate random number generator, necessary to 
implement one of the three methods described in this study. 
Miye third portion consists of the standard Gaussian random 
member generator and those steps necessary for generating 
Samples Of the random vector, X, such that a matrix is 
meeled with one thousand random vectors. Two times were 
considered when evaluating the programs: the "set-up" time 
and the "repetition" time. These times correspond to 
SeecUtion time for the second and the third portion respec- 
tively of the computer program. Similarly, the memory 
Space requirements shown in the tables make use of Two 
numbers. The first number provides an indication of the 
mourney Ol Space required for the computer programs for 

each of the three methods. This number includes the sub- 


routine EIGEN in the rotation method and the function GRN 
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in all cases. The second number is the total space 

required for the program plus any external function used 
such as square roots, absolute values, exponentials, input 
and output devices, etc. The second number varies as to 

the computer system being used and in this case the computer 
was the IBM 360/67. The headings on two columns of Tables TI, 
IT, and III are labeled "Specific Conditioning" and 
"Specific Triangularization" respectively. These refer to 
Pomputer programs written specifically for 2 x 2 and 3 x 3 
covariance matrices, and are adaptations of the general 
programs. 

Bescon uMemcdaramconwarmed in tables ly its eand Til, 
the best method to use for generating random vectors from 
covariance matrices dimensioned greater than 3 x 3 appears 
wo be the triangularization method. This method requires 
less memory space and considerably less execution time than 
the other two methods. If the user is interested in 
generating samples from a bivariate or trivariate normal 
distribution, then either a specific conditioning or triang- 
ularization method can be used. Both methods require 
about the same amount of execution time and their space 
requirements are essentially the same. 

All methods use, as a basis, the univariate normal 
random number generator. Therefore, in order to establish 
a lower bound on times required in generating normal 
random vectors, the time required to generate normal random 


numbers was obtained. These numbers should serve as lower 


1ef/ 


bounds on the times required to generate a sample of one 

thousand repetitions of a 10, 5, 3, and 2 element vector 

respectively. The lower bounds obtained in this way are: 
10,000 numbers - 1.06740 seconds 

0.53209 seconds 


0.32302 seconds 
0.20332 seconds 


5,000 numbers 


3,000 numbers 


2,000 numbers 


A specialized covariance matrix, (identity matrix), 
was used as an input to all three methods to provide the 
least cumbersome, mathematically speaking, of all matrices 
meme wOoula be Encountered. Ihe results of this test are 
tabulated in Table III. Each method maintained its overall 
ranking with respect to memory space and time requirements. 
The rotation method displayed a sharp decrease in set-up 
time which is reasonable as eigen values are easier (faster) 
to compute in this case. In general, each of the three 
methods required essentially the same amount of time for 
1000 repetitions using a typical covariance matrix as for 


the specialized covariance matrix. 
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ThyVev iia) sae 


MEMORY SPACE REQUIREMENTS 


The first number is the amount of core space required 
for the computer program and the second number is the 
total amount of core space required in execution of the 
program. This second number includes, for example, 
eAver@al funetions amd input—-ouvous devices. Each 


Space requirement is in bytes. 


Matrix aa : Triangular- : 


Hs, 952743 0.7 25 


16 ,008/40, 344 


24 192/48 ,528 


4H 5992/69, 328 





IY, 


TABLE II 


EXECUTION TIME (TYPICAL COVARIANCE MATRIX) 


The first number is the set-up time, and the second 


number is the repetition time. Each execution time is 
ie Seconds. 
Matrix 


Ze KX. 2 .00273/.339851 .00106/.33800] .00496/. 37666 


.00145/.52411 


.002477.99554 







DOOWEy/ 255305 .01398/.62023 


Te) 
~ 
UO 





O3CN57 1.06255 Oey 7 lave Seem 


WN 
» 
U1 


= 
[<) 
ms 

a) 
2, 








.00795/2.52280 |. 41652/3. 72289 





70306/2.75560 


Soe clmlCc  COMd a Iopect tic Trt ange 


i Ae 
eK oe -00091/.25246 | .00101/.27487 font 


Se ae 2OOMO1/.39530) .00M097 .39572 


: 
| 


AG edinoe Shiba 


BABCUITON Time (SPECIALIZED COVARIANCE MATRIX TI) 


the first number is the set-up time, and the second 


number is the repetition time. Each execution time is 
in seconds. 


weds Conaa tona ne ESI eu Rotation 
S1Ze Za On 


22x 2 .00278/7. 34299 | .00122/. 33509 | .00171/. 36376 













Box. S OO O27 sebood | s.O0lGw/, 512231002347 160250 
px 5 pega t4 706025) 2002607 93449 10040371, 19545 
mao 57260933 ( 0073372. 426451 01162/3. 416i 


—— ao a ieee 
.00104/7.25215 | .00078/. 24921 
.00135/7. 37903 | .00096/. 38795 ss 
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APPENDIX A 


PISCUSSEON OF THE COMPUTER PROGRAMS 


All three computer programs were written in Fortran IV. 
They are not optimal, i.e., improvements to the programs 
MeaemmoOoss Lble. thew results ol these improvements could 
change the value of memory space requirements and program 
execution times, but should not change the comparative 
ranking among the three methods. 

Two univariate normal random number generators that 
were available for use in the present study were GAUSS and 
GRN. The subroutine GAUSS is part of the IBM scientific 
subroutine package and is based on the mean value theorem. 
The function GRN was compiled by the Naval Postgraduate 
eemool staff, and is based’ on the G. Marsaglia technique 
[Ref. 7]. The execution time of GRN is approximately ten 
emis faster than GAUSS, and for this reason it was used 
in all the programs. The rotation method uses the sub- 
routine EIGEN which is also available as part of the IBM 
scientific subroutine package. Three subroutines, (DTERM, 
REDMAT, and COFACT), were developed for use in the con- 
fer onings method. The subroutine DTERM computes the deter- 
minant of an n-square matrix. The subroutine COFACT 


,th -th 


removes the elements from the row and j 


Column er -an 
n-square matrix (Z), and provides the (n-1)-square matrix 


thal 1S required in the computation of the cofactor Osa: 


Ze 


The subroutine REDMAT reduces a matrix to the desired size. 


(k) 


iets capable of providing the user with a 2 from “an 
fen 1) Cimension matrix where k < n. 

Each method was evaluated using four different size 
Mavrices as listed in Tables I, II, and III. These were 
selected as being typical of the range of matrix sizes that 
iene be encountered in practice. However, the programs 
were written in general terms so there are no program 
imposed restrictions on matrix size. 

As mentioned previously, specific computer programs 
were written for the 2 x 2 and the 3 x 3 covariance matrices 
using the conditioning and triangularization method. This 
was done in an attempt to evaluate the possible reduction 
in memory space and execution time under those obtained 
using the general programs with input dimensions of 2 and 
3. Of the three methods considered, these two appeared to 
Be adaptable to simpler, more concise programs for these 
small dimensions. Writing specific programs for covariance 
matrices dimensioned greater than 3 x 3 would become dif- 
ficult, and a general program would probably have to be 
used. No attempt was made to write a specific program for 
meee rotation method since computation of eigen values and 
vectors are involved, and there were no apparent means of 
meomcing the computation time under that of using the general 


Subroutine EIGEN. 
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ine computer programs, used in generating multivariate 
normal random vectors for each of the three methods, are 


mreluded aS a portion of this study. 
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COMPUTER PROGRAMS 
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